On the survival of the quantum depletion of a condensate after release from a magnetic trap

We present observations of the high momentum tail in expanding Bose–Einstein condensates of metastable Helium atoms released from a harmonic trap. The far-field density profile exhibits features that support identification of the tails of the momentum distribution as originating in the in-situ quantum depletion prior to release. Thus, we corroborate recent observations of slowly-decaying tails in the far-field beyond the thermal component. This observation is in conflict with the hydrodynamic theory, which predicts that the in-situ depletion does not survive when atoms are released from a trap. Indeed, the depleted tails even appear stronger in the far-field than expected before release, and we discuss the challenges of interpreting this in terms of the Tan contact in the trapped gas. In complement to these observations, full quantum simulations of the experiment show that, under the right conditions, the depletion can persist into the far field after expansion. Moreover, the simulations provide mechanisms for survival and for the the large-momentum tails to appear stronger after expansion due to an acceleration of the depleted atoms by the mean-field potential. However, while in qualitative agreement, the final depletion observed in the experiment is much larger than in the simulation.

In early tests of our measurement sequence we noticed a contamination of the signal by spurious counts. Specifically, an extremely dilute remnant of the m J = +1 cloud was observed as a small peak in counts in the region −10µm −1 k z −8.5µm −1 . We inferred these were remnant counts from the m J = +1 cloud as they were still visible when we ran an experimental sequence without the Landau-Zener transfer, but including the Stern-Gerlach pulse.
While the cause of the appearance of the remnant counts is unclear, we observe that the count density outside the region of interest is similar in both the shots with the RF pulse (Landau-Zener transfer) and those without. We note that only about one in a million atoms from the m J = 1 cloud are present in this manner in a given shot. On average, this accounted for approximately N mJ=+1 = 1 atom per shot when running the sequence without the Landau-Zener transfer, and approximately half this amplitude after applying the pulse. As the signal of interest was not the angular density profile but the integrated number of counts, we corrected for this contribution by subtraction, N k min ,k max = (N detected k min ,k max − η 1 N m J =+1 )/(η 0 × η Q × Ω ROI ), where Ω ROI is the fraction of the sphere retained by the windowing procedure described in the main text. and η Q is the detector quantum efficiency. Such counts constitute about 10(5)% of the detection events in the ROI and are not sufficient to account for the excess tail strength. We hypothesize that the remnant counts are atoms transferred into the m J = 0 state by non-ideal behaviour of the Stern-Gerlach pulses or magnetic field switches. Figure S1. Determining the RF transfer efficiency. The time-of-flight profiles of each pulse are resolved (a) by applying a weak Stern-Gerlach pulse during the time of flight. The pulses are aligned with respect to their centre-of-mass (b) and used to determine the pointwise fraction ((c), dotted line). Detector saturation is evident in the peaks (dashed lines), but not in the thermal tails (solid lines), which are used to compute the transfer efficiency. Because of its lower flux, the m J = −1 pulse does not show any clear evidence of saturation (d) and is used to determine the thermal fraction and hence N 0 . Table S1 shows the absolute and relative tail amplitude estimates Λ fit and Λ fit /Λ pred , respectively, that proceed from variations of the ROI and assumed quantum efficiency in the analysis of the experimental data. It is seen that neither the uncertainty in the collection efficiency nor the choice of elevation angle cutoff φ c have a significant effect on the findings. QE fit r 2 Λ fit Λ fit /Λ pred 0.05 0.83 0.2(0.1,0.3) 6.8(4.5,9.1) 0.06 0.83 0.3(0.2,0.4) 7.4(4.9,9.8) 0.07 0.83 0.3(0.2,0.4) 7.8(5.2,10.5) 0.08 0.83 0.4(0.3,0.5) 8.3(5.5,11) 0.09 0.83 0.5(0.3,0.6) 8.7(5.8,11.6) 0.1 0.83 0.6(0.4,0.7) 9.0(6.0,12.1) 0.11 0.83 0.6(0.4,0.8) 9.4(6.2,12.5) φ c fit r 2 Λ fit Λ fit /Λ pred 80 • 0.82 0.1(0.0,0.1) 9.3(6.0,12.7) 70 • 0.86 0.2(0.1,0.3) 9.2(6.4,12) 60 • 0.83 0.4(0.3,0.5) 8.3(5.5,11) Table S1. Tail predictions with modified ROI or assumptions about the value of quantum efficiency (QE). A change in quantum efficiency (QE) presents by changing the value of ε used in the prediction and also through the factor of N 7/5 0 used to compute the condensate density n 0 . These effects partially cancel to produce a weak scaling of Λ fit /Λ pred with respect to ε. Re-running the analysis using different QE and the standard ROI with φ c = π/3 yields fits that barely differ in the goodness-of-fit criterion and present comparable results for Λ fit . We also find that the choice of collection area defined by the elevation angle cutoff φ c (while using the best QE value of 0.08) has a weak effect on the result, but below statistical significance. Terms in brackets are the upper and lower 95% confidence intervals. Neither the uncertainty in the collection efficiency nor the choice of elevation angle cutoff φ c have a significant effect on the findings.

Statistical details
In Figure S2 we illustrate how fixing the scaling exponent α when fitting to a given dataset obscures the enormous covariance between α and the scale coefficient C α . The red diamonds in Figure S2 show the fit exponent and amplitude obtained from a fit of Eqn. (13) to single data sets, each representing a different BEC density. The variation in the fit exponent is small (mean 4.2, standard deviation 0.4), and although the mean α from the fits is only half a standard deviation (one standard-error interval) different from 4, the corresponding difference in amplitude C varies exponentially with α, spanning over six orders of magnitude. Figure S2. Illustrating the large systematic errors in a fit to the density using a power law ansatz. If the fit to the density profile (Fig 1 in the main text) ) is replaced by Eqn. (13) (also in main) and α and C α used as fit parameters, it gives a wide variation in best-fit exponents and scale coefficients between each data set (red diamonds). The red square shows the mean fitted α and (geometric) mean C α with standard deviation (in α) shown as error bars. Blue circles show the amplitude coefficient C α obtained from the same data sets with α fixed at a number of values within the range of the free fit results. The choice of α strongly determines the coefficient C α , but the error bars (standard errors in fit parameters) are smaller than the markers in all cases.
The blue circles in Figure S2 show the amplitudes returned from fits to data from all datasets when α is constrained to the set values as shown. These set values are chosen within about one standard deviation of the mean α from the two-parameter fits to separate datasets. These blue fits scarcely differ in their goodness-of-fit criterion (the mean square error) and so offer no obvious way to reconcile the expected distribution with these divergent statistical conclusions. Furthermore, this is not reflected in the error estimates in the fitting routines: The error bars representing the uncertainty in parameters from the fixed-α fits are smaller than the markers used in Figure S2. A linear fit reveals that d log 10 C α /dα ≈ 6.8.
Ultimately, the question of whether the data is drawn from a power law at all is not amenable to a decisive conclusion. For example, a log-normal distribution can produce similarly accurate predictions to the power law (see Figure 2 in the article) although there is no physical hypothesis that predicts such a distribution. This points to one of the primary challenges with power laws; the exponents are strongly entwined with the rate of occurrence of rare events, which by definition are subject to large statistical fluctuations and and thus subvert even the most meticulous investigations. These problems with fitting power laws are ubiquitous, and made more difficult by the small range of k which are visible in the helium experiments. In general, estimating the exponent of a purported power law is difficult and requires data spanning several orders of magnitude in scale 59, 60, 84, 85 , which are not present in either of the helium experiments to date.
Nevertheless -the analysis of Fig. S2 does allow us to at least say that the data is inconsistent with power laws outside the range 3.3 < α < 4.7.

Observables
Detailed data extracted from the simulations are given in Table. S3.
The total quantum depletion of the condensate δ B is given by

S3
Stochastic averaging over all trajectories in the ensemble is denoted by · stoch. . The density of the depleted particles is evaluated by the standard positive-P expression, The density of condensate is augmented by the 1 − δ B factor when calculating observables to conserve overall particle number. The k-space fields here are normalised as Shown is the calculated contact C sim at the end of the initial state generation (t = 0), as a function of the quench ramp time t ramp . The horizontal dashed line shows the LDA prediction C = 0.102 pm −1 , the blue data point the ensemble that was deemed to agree, and was used for subsequent simulations for t > 0. Lower panel: Evolution and stabilisation of the contact. The blue shaded area shows the duration of the ramp from t start to t start + t ramp . Orange shading denotes the error bars. S = 4000 trajectories in all cases.

Initial ensemble
A procedure for generating the equilibrium state has been developed in the Wigner representation 32, 86-89 . Unfortunately we cannot use this directly, nor a direct Wigner representation of the quasiparticles because the condition required for correctness of the Wigner representation -that there are more particles than modes -is very far from being met 33, 90, 91 (In fact here we have about O(1000) particles in the depletion, and O(10 6 ) modes.). It is also unclear how to translate a local density formulation of depletion in a uniform section of gas to a positive-P ensemble without introducing discontinuities. Instead we turned to dynamically generating a state with the required quantum depletion. We begin with a fully condensed ground state with ψ B = ψ B = 0 and φ (x) = φ 0 (x), the ground state of the Gross-Pitaevskii equation (21). The latter is obtained by imaginary time propagation of the GPE augmented with an appropriately chosen chemical potential µ according to (8) which sets the central density.
Our first attempt to generate the equilibrated quantum depleted state thus started with φ 0 (x) and then adiabatically ramped the interaction from g = 0 to the experimental value while evolving the equations (21)-(23). This did not work for two reasons. Firstly, a very strong collective oscillation was induced, since the width of the Thomas-Fermi ground state depends strongly on g. Secondly, very long time evolutions succumb to excessive noise in the positive-P simulation and produce a state that is too noisy to be useful. Note that in the positive-P representation different ensembles can represent the same state but exhibit very different noisiness and practical usefulness 92, 93 .

S4
The second attempt began with the Gross-Pitaevskii ground state and the target physical interaction strength g in the GPE (21), while slowly ramping the interaction strength in the Bogoliubov equations (22)-(23). This eliminates the main oscillations in the mean-field evolution. The least unwanted nonadiabatic disturbance occurs when g is ramped only within the projected part of the Bogoliubov evolution as per with g B (t). However, over adiabatic timescales, this still introduced far too much noise in the Bogoliubov fields ψ B , ψ B to be useful.
To work around the problem, we take advantage of the fact that an instantaneous quantum quench of a weakly interacting condensate produces depletion with a similar time-integrated momentum profile dt n B (k,t) ∝ 1/k 4 to the ground state value, but with a somewhat higher depletion, as described in detail in 61,80 . Simulations show that non-instantaneous but rapid quenches produce lower amounts of quantum depletion, as shown in Figure S3 (left panel). To generate these ensembles, we used the following ramp: starting at negative t start , and evolving till t = 0. Thus t ramp was the ramp time, and the remaining time: |t start | − t ramp an equilibration time. We carried out a calibration like that shown in Figure S3 (left panel) for each set of cloud and trap parameters needed. Then for the initial in situ ensemble in the trap we chose the ensemble generated with the t ramp that produces total quantum depletion in agreement with the in situ value (9) for the condensate ground state (like the blue data point in Figure S3, left panel). The value of t start was −101µs for the weakly trapped cases (ω = 45 × 425 × 425 Hz andω = 201 Hz) and −198µs for the strongly trapped cases (ω = 71 × 902 × 895 Hz andω = 393 Hz) detailed in Table. S3. Figure S3 (right panel) shows an example of the evolution of the contact in situ during this initial state generation. Properties of the initial states are labelled (CT) and (ST) in Fig. 3 (a) (main text).

Determination of the contact and tail strengths
Much as in the experiment, the k-space density n(k) in the simulations is very noisy in the asymptotic region of large |k|, and a lot of averaging is needed to extract the contact. In light of the discussion on power laws, the exponent α was not fitted here, but we did calculate both the fit coefficients C sim on the function n B (k) = C sim /k 4 with α = 4 assumed and the count of particles in the tails, N k min ,k max as in Table S2, in order to have a systematic view of how correlated the two kinds of results are. We proceed as follows: The simulations provide a density of depleted particles n B (k). Like with the experimental counts, we keep only density that is far enough away (φ c , usually 60 o ) from the long axis of the cloud. The simulated system is axially symmetric around the long axis. We do not restrict counting to the vicinity of the vertical axis because here there are no detector irregularities to avoid, and this allows us to improve the signal to noise ratio.
The particle counts N B (k) = n B (k)/V are finely binned according to the absolute value of the momentum k = |k|, giving total bin count N k . The field is simulated on a square lattice in k space, in which each site corresponds to a k-space volume of ∆V K = (2π) 3 /V with real space volume V . Therefore the total bin volume V k is obtained by binning the site volumes, and the mean density in each bin is N k /V k so n k = (2π) 3 N k /V k . Each bin gives an estimate of the corresponding apparent contact: The statistical error estimate on c(k), ∆c(k), is obtained by averaging subensembles, then using the central limit theorem (CLT) on the subensemble averages. Figure S4 shows example values and error estimates, while also demonstrating the difficulties involved.
We can see that for a significant range of k values the c(k) estimate gives fluctuations around a constant value. However, several difficulties in extracting the overall trend are also evident. Low k values are not representative because the particles in the shaded region never emerge from the condensate and are not measured in the experiment. Therefore we remove data below a value k inner from consideration. k inner is chosen such that it does not include any depletion that would significantly overlap with the BEC in the final expanded cloud and be obscured by it. This corresponds also to the energy at which the quasiparticle spectrum becomes particle-like, since the mean field energy in the centre of the trapped cloud is responsible for both effects.
High k values, on the other hand, suffer from statistical error sufficiently high as to make them useless, and so we also choose a maximum k outer value, and only use k ∈ [k inner , k outer ] to extract the contact estimate. The extent of the useless high k region changes with time during the simulation. To systematically adapt our fitting region to this, we choose k outer according to the calculated statistical error such that and ∆ max is chosen once for all times in a simulation. This gives a time-dependent fitting region k ∈ [k inner , k outer (t)]. The final contact estimate C sim (t) is the mean of all c(k,t) values in the fitting range. Error bars on C sim (t), ∆C sim (t), are obtained from a CLT estimate of the error in the mean of the points used, after binning enough neighbouring points to encompass the autocorrelation (i.e. so that neighbouring bins are uncorrelated). The ∆ max is chosen low enough so that the statistical error in individual points ∆c(k,t) at the large k end does not put off the overall contact estimation C sim (t) ± ∆C sim (t). Operationally we choose the value of ∆ max at which the error introduced by adding successive high k points introduces more uncertainty than the reduction thanks to a larger ensemble. We also found that the relative factors between Eqn. (14) and the simulated tails (i.e. C sim /C and A sim ) depend on choice of cutoff angle φ c . We find that apparent tail strengths C sim and A sim are larger for smaller collection regions that are more tightly concentrated about the strong trapping axis, whereas larger collection angles (that include areas closer to the weak axis) produce lower apparent tail strengths. Data on A sim are shown in Table S2 while C sim was 1.6(1), 1.9(2), and 2.2(3) times C for φ c values of 60, 70, and 80 degrees, respectively. n o N k min k max φ c N k min ,kmax ratio A sim (10 6 µm −3 ) (µm −1 ) (deg) final final/in situ Rapid release from trap (CE) 5.237(16) 2.0 3.5 60 193 (8) 1.38(5) 12.09(6) 2.25 3.5 60 362(10) 1.43(4) 6.86 (3) Table S2. Tail strength data in simulations. Based on final times in the simulations described in Table S3, referenced by the value of n 0 N. N k min ,k max is calculated as per Eqn (14) in the main text. The ratio A sim of tail strength in the expanded cloud compared to Tan theory predictions is obtained by dividing N k min ,k max at the final time by its value in situ at t = 0. k max is chosen to still contain all 4π steradians inside the square lattice, while k min to avoid overlap with the expanding condensate.  Table S3. Main simulation data and parameters for fitting of k −4 tail amplitude C sim . The time t for which data are calculated is counted relative to the start of the trap release. Abbreviations (CT,CE,. . . ) as in Figure 3a of the main text. The range k inner -k outer here was chosen as explained in the text and used to obtain the estimate and uncertainty for C sim . In all cases, S = 4000 trajectories averaged. N B is the number of Bogoliubov field particles as per (S1).

S7
Toy two-mode model of escape In a uniform gas in the Bogoliubov approximation, the a k mode is coupled only to a † −k and the condensate. The Bogoliubov-de Gennes equations for these modes can then be written 81 d dt Notice the provision for a time-dependent background condensate density n(t). We stay in the real-particle basis rather than the Bogoliubov quasiparticles b k which allows to avoid calculating time-dependent changes of the coefficients u k , v k . The chemical potential is µ(t) = gn(t) 81 . The equations (S9) can be used as input for equations of motion of the low order moments ρ(k) = a † ±k a ±k , A(k) = a k a −k , and A * (k) = a † k a † −k , for example dρ(k)/dt = (d a † k /dt) a k + a † k (d a k /dt) . Assuming equal initial occupations ρ(k, 0) = ρ(−k, 0), one obtains an evolution equation for two coupled quantities ρ(k,t) and A(k,t) = A r (k) + iA i (k). It is The initial conditions are ρ(k, 0) = v 2 k , A(k, 0) = v k u k , corresponding to the Bogoliubov ground state. Taking n(t) constant one obtains the solution (16). A more realistic model is obtained by estimating the n(t) that particles in the k mode experience as the trap is released and the condensate expands. Hydrodynamic expansion of a condensate with a Thomas-Fermi profile suddenly released from the trap 94 gives a self-similar expansion with widths approximated by W j (t) = W j (0) 1 + (ω j t) 2 . The initial widths being the Thomas Fermi radii in the trap R j = (1/ω j ) 2gn 0 )/m, j = {x, y, z}. From this, and assuming ω y = ω z = ω ⊥ =ωλ 1/3 (approximately true in our setup) the condensate density at any point in the cloud decays as n c (r,t) = n c (r, 0) Escape is much easier, however, for atoms near the edge of the cloud than in the centre, because they can avoid reabsorption by faster free-flight from the high density region. We categorise the initial position of the studied section of the gas with the quantity R 0 = y 2 + z 2 R ⊥ = ω ⊥ y 2 + z 2 2gn 0 /m .
The initial density is then The motion of the atoms is estimated by assuming they travel in the outward direction at velocityhk/m relative to the hydrodynamic expansion (which also accelerates them to some degree). Hence their position is estimated as The condensate at this position comes from hydrodynamic expansion of the initial condensate at position r c0 (t) = r(t) (S15) The local density at position r at time t is then n c (r c0 (t),t) according to (S11). Finally, to roughly include ramp time in the analysis, consider the density in a trapped gas that adiabatically follows the trap evolution ω j (t) = ω j (0)e −t/τ release . The overall density for a piece of gas that stays in the same relative position in the cloud will then follow n adiabatic c (t) = n(0)e −3t/τ release . (S16)

S8
To incorporate this effect into the model it is combined (somewhat ad hoc) with the expression n c (r c0 (t),t) for instantaneous release to give the final estimate of local density for the k momentum escaping atoms used in (S10) as n(t) = n(0)e −3t/τ release + (1 − e −3t/τ release ) n c (r c0 (t),t). (S17) with the help of (S11) and (S15).